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Abstract 

Impact craters are among the most studied geomorphic planetary features because they yield information about the past geological 
processes and provide a tool for measuring relative ages of observed geologic formations. Surveying impact craters is an important 
task which traditionally has being achieved by means of visual inspection of images. The shear number of smaller crater present 
in high resolution images makes visual counting of such craters impractical. In this paper we present a method that brings together 
a novel, efficient crater identification algorithm with a data processing pipeline; together they enable a fully automatic detection of 
sub-km craters in large panchromatic images. The technical details of the method are described and its performance is evaluated 
using a large, 12.5 m/pixel image centered on the Nanedi Valles on Mars. The detection percentage of the method is ~70%. The 
system detects over 35,000 craters in this image; average crater density is 0.5 craters/km 2 , but localized spots of much higher crater 
density are present. The method is designed to produce “million craters” global catalogs of sub-km craters on Mars and other 
planets wherever high resolution images are available. Such catalogs could be utilized for deriving high spatial resolution and high 
temporal precision stratigraphy on regional or even planetary scale. 
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1. Introduction 

Impact craters are structures formed by collisions of mete¬ 
oroids with the planetary surface. The importance of impact 
craters stems from the wealth of information that detailed anal¬ 
ysis of their distributions and morphology can bring forth. For 
example, statistics of crater sizes form the basis for geological 
stratigraphy of the Moon and Mars [8, 29, 23], Craters counts 
are the only available tool for measuring remotely the relative 
ages of geologic formations on planets. In addition, knowl¬ 
edge of crater morphologies enables studies of a number of 
outstanding issues in planetary geomorphology, such as: the na¬ 
ture of degradational processes [20], regional variations in ge¬ 
ologic material [6], and distribution of subsurface volatiles [7]. 
Thus, surveying impact craters is an important task in planetary 
research. Presently, all such surveys are done manually. For 
example, the most comprehensive catalogs of craters on Mars 
contain information on 42,283 craters [4], and 57,633 craters 
[16], respectively. They contain only craters larger than 5 km 
in diameter. For the Moon the most comprehensive catalog [1] 
lists 8,497 named craters that are larger than few kilometers in 
diameter. For Mercury and the icy satellites of outer planets 
much less extensive catalogs of large craters have been com¬ 
piled. No comprehensive catalogs of smaller, sub-km craters 
exist for any planetary body. This is because building such a 
data set is a very laborious process, ill-suited for visual detec¬ 
tion but well-suited for an automated technique. 

Automating the process of crater detection is the only practi¬ 
cal solution for comprehensive surveys of smaller craters. Such 
surveys can deliver regional or even global coverage making 
million craters catalogs a reality and facilitating an objective 


and repeatable analysis, well beyond and above of what can be 
analyzed using manually processed data sets. 

Auto-detection of craters in planetary images is not a new 
idea; Salamuniccar and Loncaric [17] have tabulated 73 publi¬ 
cations devoted to various techniques of crater detection. We 
refer to their tabulation for a complete record of previous work 
on crater detection algorithms. Despite this large body of work, 
practitioners of planetary science continue to count craters man¬ 
ually, and we are not aware of any applications of crater detec¬ 
tion algorithms to actual planetary research. This is because the 
early approaches to crater detection were not developed beyond 
a demonstration stage and were not demonstrated to be robust 
to changes in a character of planetary surface thus limiting their 
appeal in actual applications. Only the most recent crater detec¬ 
tion algorithms, that borrow techniques from the field of face 
detection, appear to be robust enough for practical applications. 
Kim et al. [12] use combination of edge detection, template 
matching, and neural network-based scheme for recognition of 
false positives to construct a system for auto-detecting craters 
on Mars. In a different approach, Martines et al. [13] adopted 
a boosting algorithm, originally developed by Viola and Jones 
[27] in the context of face detection, for identification of craters 
on Mars. Although these new approaches show high detection 
accuracy, the underlying technology is complicated and its ro¬ 
bustness to different terrain types and to image quality were not 
demonstrated. 

In this paper we present a different approach to auto¬ 
detection of craters in panchromatic planetary images. The key 
insight behind our method is an observation that a crater can be 
recognized in an image as a pair of crescent-like highlight and 
shadow regions. Therefore, our method utilizes mathematical 
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Figure 1: Diagram illustrating consecutive steps in mathematical morphology-based algorithm for identification of craters. 


morphology [14, 18] for the task of crater detection; scale and 
rotation-invariant shape filters are used to identify crescent-like 
regions in an image. Subsequently, we utilize supervised ma¬ 
chine learning techniques to distinguish between regions that 
are craters and those which are not. This approach, based on a 
relatively simple principle, promises to be robust and efficient, 
making it well-suited for large-scale, regional or global surveys 
of small craters. Our focus is on surveys of sub-km craters be¬ 
cause larger craters can be identified either manually [4, 16], or 
using a topography-based algorithm [21, 22]. The present paper 
describes the core crater detection algorithm based on principle 
of mathematical morphology, gives the assessment of its perfor¬ 
mance, and shows how it can be utilized for large scale surveys 
of craters. 

The rest of the paper is organized as follows. Section 2 de¬ 
scribes the core crater identification algorithm. In Section 3 
we present the results of applying our method to a very large, 
12.5 meters/pixel panchromatic image of Mars taken by the 
High Resolution Stereo Camera (HRSC) instrument on-board 
the Mars Express spacecraft; both accuracy assessment and full 
image crater survey are presented. Conclusions and future re¬ 
search directions are given in Section 4. 

2. Methods 

Fig. 1 shows a diagram illustrating individual steps in our 
core algorithm. The input is a panchromatic image which con¬ 
tains many highlight and shadow features. The algorithm pro¬ 
cesses highlight and shadow features in parallel using inverted 
image to process the shadow features. The goal of the chain 
of processes is to eliminate all features that are not indicative 
of craters. The remaining highlight and shadow features are 
matched to each other to flag regions that are candidates for 
being craters. In the final step a supervised machine learning 


technique is used to distinguish between true craters and non¬ 
crater regions. The feature processing pipeline uses extensively 
the attribute filters [5] including power filter, area filter, and 
shape filter. An attribute filter uses a criterion evaluated us¬ 
ing attributes of image features to decide whether to remove or 
keep unchanged any given feature. For example, if an attribute 
is the number of pixels contained in a feature and the criterion 
is “attribute value > A”, then a corresponding attribute filter re¬ 
moves all image features that contain less than A pixels. The 
attributes used in our method are invariant to translation and ro¬ 
tation, which means that corresponding attribute filters remove 
all image features that do not satisfy the criterion regardless 
of the position and orientation of each feature. The shape at¬ 
tributes used are also scale invariant. Thus, a single application 
of an attribute filter to an image removes all undesired features 
thereby saving computational cost in the processing pipeline. 

2.1. Preprocessing 

An image of planetary surface contains many highlight and 
shadow features that are either not parts of craters or are parts 
of large craters that are of no interest to our present method. 
Thus, in the preprocessing step, we apply a series of filters to 
remove such features. 

• Background removal: This step removes background fea¬ 
tures such as mountains that are too large to be part of the 
craters. A median filter is applied to the input image, re¬ 
sulting in an image containing only the global features of 
the background; by subtracting this median-filtered image 
from the input image, a new image is obtained which no 
longer contains the large background features. A median 
filter used in test sites uses a 201 pixels wide circular win¬ 
dow. 

• Power filter: The power filter [30] removes features that 
are (almost) invisible: they are too small and/or they lack 
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sufficient contrast to be noticed by humans and can there¬ 
fore also not be used for reliable processing. We imple¬ 
ment the power filter as an attribute filter with the power 
attribute defined as P — A(h a - hi,) 2 using variables com¬ 
puted from an image feature: its area A, its gray-level h a 
and the gray-level hj, of the darkest neighboring feature 
that is still brighter than h a . Features with / J < Pq are 
removed. We use Po = 1000 in test calculations. 

• Area filter: The area filter is an attribute filter using the 
number of pixels in a feature as its attribute. This filter 
removes all features that are considered too small for reli¬ 
able crater detection. Features with A < Aq are removed. 
We use Ao = 30 in test calculations. 

2.2. Shape filters 

After removal of features that are either too large, too small, 
or lack sufficient contrast, the next step in our method is to re¬ 
move features that have shapes irreconcilable with craters. For 
this purpose we employ a shape filter that uses a shape attribute. 
In order to achieve a sufficient level of discrimination between 
crater and non-crater features, the shape attribute is a vector 
containing Hu’s seven moment invariants [9], These invariants 
are easy to implement and efficient to compute [24], So-defined 
shape attribute is invariant to translation, rotation, and scaling 
[25]; the corresponding shape filter preserves or removes im¬ 
age features based exclusively on their shapes without regard to 
their size, orientation or location. Thus, only one pass over all 
features is needed to remove all the features whose shape does 
not satisfy the given criterion. 

Shape rejection criterion is based on the minimum Euclidean 
distance between a shape attribute of a feature and a set of 
shape attributes belonging to a library of reference shapes cho¬ 
sen to represent characteristic highlight and shadow regions of 
the craters. In test calculations the number of reference shapes 
for the highlights and shadow regions was 17 and 12, respec¬ 
tively. An image feature having distance of less than 0.05 to at 
least one of the reference shapes is preserved for further pro¬ 
cessing. 

2.3. Matching highlight and shadow regions 

Preprocessing and application of shape filters yield two im¬ 
ages: one containing highlight regions that are likely to be the 
parts of craters and another containing shadow regions that are 
likely to be the complementing parts of craters. The next step in 
our method is to match highlight regions to shadow regions so 
each pair corresponds to a single crater. This step is important 
because it incorporates additional knowledge into final crater 
identification. As each crater must appear as a pair of highlight 
and shadow regions with proper sizes, shapes, and relative lo¬ 
cations, regions that cannot be matched are not considered to 
correspond to craters. To obtain crater candidates we match the 
highlight and shadow regions using simple rules. Regions that 
cannot be matched are discarded. Each possible pair of high¬ 
light (H) and shadow (S) regions are considered. A pair (H, S) 
must satisfy all of the following rules to be considered a crater 
candidate: 


1. Distance between H and S must be smaller than a threshold 
proportional to the regions size. This removes pairs of H 
and S regions that are too far away from each other. We 
use 1.65 sfAfi as the threshold, where Ah is the area of H. 

2. The regions H and S should have similar sizes. In present 
calculation the difference in size cannot be larger than 4. 

3. The elongation of combined H and S regions must be 
smaller than a threshold; only pairs with a round combined 
shape are preserved. The measure of elongation is pro¬ 
vided by the first of Hu’s seven moments invariants. We 
reject pairs with combined shape having the first moment 
invariant larger than 3. 

4. The elongation of combined H and S regions must be less 
than individual elongations of H and S. 

5. The H and S regions must align along the orientation dic¬ 
tated by a solar illumination angle. 

2.4. Supervised classification 

The above rules have been established experimentally and 
they provide a “first cut” list of crater candidates. For further 
refinement of that list we use supervised classification. First, 
the shape and size of each candidate is estimated by filling the 
space between its constituent H and S regions using a morpho¬ 
logical closing operation with an ellipse as the structuring el¬ 
ement. Each candidate is assigned a multidimensional vector 
of attributes which consists of Hu’s seven moments invariants 
computed for the approximate shape of the candidate and the 
variables computed in the process of executing the rules enu¬ 
merated in section 2.3. We discriminate between craters and 
non-craters among the candidates using the technique of super¬ 
vised classification. A supervised learning uses a training set 
of crater candidates, for which crater or non-crater labels have 
been already assigned manually by an analyst. This training set 
is utilized by a learning algorithm to construct a labeling func¬ 
tion (a classifier) that is then applied to label all other crater can¬ 
didates. Such classifier is an extensive, computer-derived set of 
rules that reflects a connection between candidates numerical 
features and their labels. In this paper a classifier is constructed 
using the C4.5 decision tree technique [15], Applying the re¬ 
sultant classifier to the list of crater candidates yield the final 
catalog of identified craters. 

3. Results 

3.1. Test site 

We have selected the HRSC nadir panchromatic image 
h0905_0000 to serve as a test site. This large (8248 X 
65448 pixels) image represents a significant challenge to au¬ 
tomatic crater detection because of its size and uneven illumi¬ 
nation. The site is located between longitudes of -47°.22E and 
-49°.0E, and the latitudes of 0°.3S and -18°.48N. Fig. 2 shows 
the geographical and geologic context of the test site; the site is 
located in the Xanthe Terra, centered on Nanedi Valles and cov¬ 
ers mostly Noachian terrain. With image resolution being 12.5 
m/pixel, site’s north-south extension is 818 km and its east-west 
extension is 103 km. 


3 




Figure 2: Geographical and geologic context of the test site. The HRSC image 
h0905_0000 is overlaid on a map showing regional topography and geological 
units. The location of the training tile is indicated by a red rectangle, and the 
location of the test tile is indicated by an orange rectangle. 


A file containing the h0905_0000 image has a size of 515 
Mb. Because of memory considerations, a crater detection al¬ 
gorithm cannot process the entire image at once while running 
on a typical PC computer. We have developed a processing 
pipeline for surveying craters in high resolution images (like 
the h0905_0000 image) having very large file size. First, we 
tessellate the entire image with overlapping tiles, so each tile, 
having relatively small size, can be processed separately. The 
tiles overlap to ensure detection of craters located on the bound¬ 
aries of the tiles. The h0905_0000 image was tessellated into 44 
x 6 = 264 tiles, each containing 1700 x 1700 pixels. Neighbor¬ 
ing tiles have a 200 pixel overlap. Craters are detected from 
each tile and a local catalog of craters identified within a tile 
is created. By combining all tile-specific catalogs, a complete 
catalog for the entire image is created. Finally, duplicate craters 
resulting from tiles overlapping are eliminated. 

3.2. Training and testing tiles 

Before applying our crater detection algorithm to the entire 
h0905_0000 image we first apply it to a couple of single tiles 
having dimensions of 1700 x 1700 pixels each. We have chosen 
a tile labeled 3224 as a “training site” and a tile labeled 3225 as 
a “testing site.” The tile label reflects horizontal and vertical 



Figure 3: (Top) A panchromatic, 12.5 m/pixel resolution image of the tile 3_24; 
the red arrow points toward a 1 km size crater to give a sense of scale. (Bottom) 
Crater candidate shapes are shown in blue overlaid on an image of tile 3_24. 


locations of the tile in a tessellated image. The location of the 
training tile is shown by the red rectangle on Fig. 2, and the 
location of the testing site is shown by the orange rectangle. The 
training site is used to manually label characteristic examples of 
crater and non-crater image features. The testing site is used to 
evaluate the performance of a classifier on an image that is not 
contributing any features to the training set. 

The upper part of Fig. 3 shows an image of the 3224 tile. The 
red arrow indicates a 1-km-size crater to give a sense of scale. 
Great majority of craters within the tile are sub-km craters. The 
lower part of Fig. 3 shows 1543 crater candidates identified by 
our algorithm using shape filters, but before applying a super¬ 
vised classification. It is clear that shape filters alone are not 
sufficient to identify the craters; they allow many false positive 
detections. 

In order to apply a supervised classification we need to es¬ 
tablish a training set, a collection of crater candidate shapes that 
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Figure 4: (Top) Training shapes: yellow outlines indicate craters, blue outlines 
indicate non-crater feature. (Bottom) Craters detected by our algorithm are in¬ 
dicated by red outlines. Note that the largest crater in the image is not identified 
because it exceeds the 5 km size limit. 


Figure 5: Comparison of craters detected by our algorithm (red) and manually 
marked by an analyst (yellow). (Top) Training site, tile 3_24. (Bottom) Test 
site, tile 3_25. 


are labeled by an expert to be either craters or non-craters. From 
among the 1543 crater candidates identified in the 3_24 tile, we 
have marked selected 69 shapes that are indeed craters. These 
positives examples are shown as yellow outlines on Fig. 4. We 
have also marked 59 features that are in the group of crater can¬ 
didates but are not craters; these are shown as blue outlines on 
Fig. 4. The positives and negatives examples are used to train 
the C4.5 classifier. The resultant classifier was applied to all 
1543 candidates resulting in marking 185 of them as craters 
and the rest as non-craters. The identified craters are shown as 
red outlines in Fig. 4. 

In the test site the algorithm finds 2112 crater candidates. The 
supervised classifier (trained on the training site) labels 234 of 
them as craters and the rest as non-crater features. 


3.3. Testing performance of crater auto-detection algorithm 

The performance of our algorithm is evaluated using training 
and testing tiles. Craters in both sites were manually marked to 
provide ground truth used for quality assessment. The ground 
truth data contains 351 craters in the training site and 360 
craters in the test site. Fig. 5 shows craters marked by an analyst 
and craters detected by our algorithm in training and test sites, 
respectively. Table 1 lists the crater counts in the both sites. The 
table is divided into two categories: all craters and the craters 
larger than 200 meters in diameter. Because our method is de¬ 
signed to detect craters larger than ~20 pixels in diameter, we 
expect only craters larger than ~200 meters to be systematically 
detected in HRSC images. However, an analyst has marked all 
craters visible to a human eye, many of them smaller than 200 
m in diameter. The fair comparison is between algorithmically 
detected craters and manually detected craters larger than 200 
m in diameter. 
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TP 

FP 

FN 

Training site (all) 

175 

10 

176 

Training site (D > 200m) 

109 

8 

45 

Test site (all) 

198 

36 

162 

Test site (D > 200m) 

120 

35 

57 



D 

B 

Q 

Training site (all) 

49.9% 

0.06 

48.5% 

Training site (D > 200m) 

70.8% 

0.09 

66.5% 

Test site (all) 

55% 

0.18 

50% 

Test site (D > 200m) 

67.8% 

0.29 

56.6% 


Table 1: Crater counts. 


Table 2: Quality of crater detection. 


We use quality factors developed in [19]: detection percent¬ 
age D = 100r P/(TP + FN ), branching factor B = FP/TP, and 
quality percentage Q = 100 TP/(TP + FP + FN), to assess the 
performance of our algorithm. Here, TP stands for the num¬ 
ber of true positive detections (identified craters that are actual 
craters), FP stands for the number of false positive detections 
(incorrectly identified craters that are not actual craters), and 
FN stands for the number of false negative detections (non¬ 
detection of real craters). D can be treated as a measure of 
crater-detection performance, B can be treated as a measure 
of delineation performance, and Q can be treated as an over¬ 
all measure of algorithm performance. 

Table 2 lists the quality factors resulting from comparing 
craters detected by our algorithm to the manually identified 
craters. This table is again divided into two categories, all 
crater and the craters larger than 200 m in diameter. As ex¬ 
pected the results are better if only craters with diameters > 
200 m are taken into account. This is because smaller craters 
lack sufficient number of pixels for reliable machine identifica¬ 
tion. Interestingly, the performance on both sites is about the 
same indicating that the limited set of examples established in 
the training site are equally representative for both sites. 

Overall, the performance of our algorithm is D ~ 70% and 
Q ~ 55 - 65%. The most direct performance comparison can 
be made with an algorithm developed by [12] which was also 
tested using HRSC panchromatic, 12.5 m/pixel resolution im¬ 
ages. They reported D ~ 70% and Q ~ 62% when tested on 
small craters in the Olympus Mons part of the h()37 OOOO im¬ 
age. Performance of D ~ 86% and Q ~ 73% was reported 
in [3] using an algorithm based on template matching. How¬ 
ever, they tested their algorithm only on images taken by the 
Mars Orbiter Camera (MOC) in its wide-angle mode and hav¬ 
ing resolution of 200-300 m/pixel. In addition, they applied 
their algorithm exclusively to sites located in the Hesperian ter¬ 
rain. In [13] the same images were used to assess a performance 
of a crater detection algorithm based on a boosting approach 
and even slightly better performance was reported. The high 
performance reported by [3, 13] may be attributed to the fact 
that large craters in a Hesperian terrain are very regular and rel¬ 
atively easy to identify. It is unknown how these algorithms 
would perform when applied for detection of small craters in 
more challenging surfaces. 

3.4. Surveying sub-km craters in the entire h0905J0000 image 

For the survey of craters in the entire h0905_0000 image 
we have used expanded training set constructed using the tiles: 
1_23, 2225, 3223, 3224, and 4225. The total number of training 
examples is 2611, with 1994 examples of characteristic crater 


shapes and 617 examples of characteristic non-crater shapes. 
We have embedded our crater identification algorithm into a 
data processing pipeline; without this pipeline it would take a 
significant time and effort to process each of the 264 tiles sep¬ 
arately. The processing pipeline employs scripts that automate 
the entire process. The supervised classifier is run using Java- 
based routines, components of the WEKA [28] environment. 
The pipeline produces the final catalog of craters and also a set 
of ArcGIS [2] project files, each pertaining to a single tile. An 
ArcGIS project contains all the information about the single tile 
in the GIS format including the tile image, crater candidates, de¬ 
tected craters, etc. These project files are valuable byproducts 
of a survey. For example, in some studies it may be crucial to 
count craters as accurately as possible in a relatively small re¬ 
gion. In such a case, an ArcGIS project file covering this region 
can be used to manually correct an automatic survey. On the 
other hand, an auto-survey offers a “carpet coverage” of sub- 
km craters over the large regions of the planet. The accuracy of 
such a survey may not be on par with the accuracy of manual 
mapping, but it is sufficient for statistical purposes and much 
more practical. 

Fig. 6 shows locations of the 35,495 craters identified by our 
process in the h0905_0000 image. For such a massive survey a 
validation process as described in the previous subsection is not 
feasible. However, visual inspection of an image with the loca¬ 
tions of the craters superimposed on it indicates performance 
on par with that reported for tiles 3224 and 3225. 

The craters identified by the auto-survey span the sizes from 
112 m to 6 km. The size-frequency distribution (SFD) of all 
the craters identified in an image is shown on Fig. 7. It is clear 
from the observed fall off in the number of craters having diam¬ 
eters smaller than ~200 m that our algorithm has a detectability 
threshold of ~200 m or about 10-20 pixels. There is also a 
smaller fall off at the upper threshold of detectability. This up¬ 
per threshold is artificially introduced to reduce computational 
time of the survey and can be adjusted. The comparison of the 
SFD constructed from detected craters with isochrons proposed 
by [10, 11] indicates the age (Noachian) of the underlying sur¬ 
face in agreement with the age inferred from counting larger 
craters [23], 

The most striking result coming out of the survey is an un¬ 
even, patchy character of spatial distribution of craters. Fig. 6 
shows a map of the local crater density as calculated using a 
moving window technique with the size of 100 pixels. This map 
shows that although an average crater density is 0.5 craters/km 2 , 
local “hot spots” with twice that density are present. The 
largest crater-enhancement spot is located just to the west of a 
a large, 43-km-size, crater centered at -47°.86E and 10°.49N. 
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Figure 6: Auto-survey of sub-km craters in the HRSC image h0905_0000. (Top) Crater centers are indicated by red dots overlaid on the image. (Bottom) Local 
crater density. In both panels the north direction is to the right. 



crater diameter (meters) 


Figure 7: Size-frequency distribution for craters identified in the HRSC image 
h0905_0000. The blue marks indicate counts using all craters in an image, the 
red marks indicate counts using only craters located within the crater-enhanced 
spots, and the green marks indicate counts using only craters located outside 
the crater-enhanced spots. Dotted lines represent isochrones, the upper solid 
line marks saturation equilibrium, and the other solid lines mark divisions of 
Amazonian, Hesperian, and Noachian eras. 


The crater density in the middle of this spot approaches 2 
craters/km 2 . The presence of these spots is not an artifact of 
the detection method. Visual inspection confirms existence of 
many more real craters in these regions. 

In order to study further these spots, we have constructed 
an inclusive mask containing these parts of the image where 
crater density exceeds 0.8 craters/km 2 . On Fig. 6 these are areas 
shown in red or darker colors. The masked region constitutes 
only 12% of site’s area, but it contains 11,362 or 32% of all 
the craters. The SFD of the craters located within the masked 
area is also shown on Fig. 7. The distribution is steeper than 
that constructed from all the craters in the site indicating rela¬ 
tively more small craters. We have also derived the SFD for the 
24,133 craters located outside the masked region. The SFDs 
constructed for these craters is not noticeably different from the 


SFD constructed for all the craters. Thus, within the spots small 
craters are even more dominant than is expected. Explaining the 
origin of these spots is beyond the scope of the present paper, 
these may be regions where secondary craters are concentrated, 
but no clear pattern of the spots is observed with the constraints 
of the image h0905_0000. 

4. Conclusions 

In this paper we have presented a novel approach to detecting 
craters in planetary images. This approach minimizes computa¬ 
tion time while achieving accuracy comparable to other recently 
published crater detection algorithms. On the h0905_0000 foot¬ 
print used (8248 x 65448 = 539,815,104 pixels) 14 hours of 
computing time were needed, which can be divided into 7 hours 
for tiling, preprocessing, and the shapefilters, 6.3 hours for the 
supervised machine learning, and 0.6 hours for removing dou¬ 
bles. The total computing time equals to an average process¬ 
ing speed of 10,710 pixels per second or 93 microseconds per 
pixel. Because of its features we propose that our algorithm 
is a viable tool for comprehensive, large-scale auto-surveys of 
sub-km craters in high resolution planetary images. The pa¬ 
per brings together a new core crater detection algorithm based 
on mathematical morphology methods with a data processing 
pipeline that allows fully automatic survey of craters in a very 
large image. 

At present an accuracy of our algorithm, as measured by the 
detection percentage, is about 70% for sub-km craters on the 
Noachian terrain. Although we did not test the performance of 
our algorithm on larger craters and/or other types of surfaces, it 
is expected that small craters on Noachian surface presents the 
toughest challenge to auto-detection. 

Future versions of our algorithm will show an improved per¬ 
formance. This can be achieved by a number of revisions. First, 
a decision to reject or retain a shape (see section 2.2) could be 
based on a machine learning instead of comparison to a limited 
number (29) of arbitrarily chosen shapes. Second, the match¬ 
ing criteria (see section 2.3) could be reviewed and potentially 
replaced completely by a decision function generated by a su¬ 
pervised classifier that uses a large training set. This means 
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combining steps described in sections 2.3 and 2.4. Finally, the 
present algorithm uses supervised classifier based on the C4.5 
method; other classification methods , such as, for example, the 
Support Vector Machine [26] should be tested and may offer 
improved performance. 

Our method can be used to expedite crater counting in order 
to establish relative and absolute chronologies of planetary sur¬ 
faces. Typically, craters with diameters > 1 km are used for 
such purpose because manual counting becomes exceedingly 
tedious with decreasing crater size. However, only relatively 
large regions can be dated by counting craters having size > 1 
km because sufficient number of craters must be accumulated 
to reduce uncertainties associated with statistical comparison 
of counts. In order to resolve various ages of smaller regions, 
counts of sub-km craters, which are orders of magnitude more 
frequent, are needed. Because counts of sub-km craters have 
lower statistical uncertainties, they can be used to recognize 
smaller differences in age, or smaller areas of different age. 
Our auto-detection method, that delivers means of compiling 
a “million craters” catalog of sub-km craters, can, in principle, 
be utilized to produce high resolution and high precision global 
stratigraphy. 

Spatially extensive surveys of small craters may produce sur¬ 
prises and lead to new discoveries. For example, our survey 
of the HRSC image h()905 0000 reveals patchy distribution of 
craters and existence of prominent spots of enhanced crater den¬ 
sity in an area that otherwise appear to be uniform in age and ge¬ 
ology. Surveying even larger regions of Mars may reveal other 
patterns in spatial distribution of sub-km craters. 

Finally, our method is pixel-based and thus independent of a 
source and resolution of an image. In this paper we used an im¬ 
age of Mars taken by the HRSC instrument, but our algorithm 
can identify craters in any planetary image. 
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